/**
 * c11_float_operations - 浮点运算：数学函数与特殊值
 * 
 * 这是《Rust基础数据类型实战Cookbook》系列的第11个案例
 * 
 * 学习目标：
 * 1. 掌握Rust中的各类数学函数及其使用方法
 * 2. 理解特殊值（NaN、无穷大、零值）的产生和处理
 * 3. 学习数学常数的使用和精度考虑
 * 4. 掌握复杂数学运算的实现技巧
 * 5. 了解数值计算中的稳定性和精度问题
 * 6. 学会处理边界情况和异常值
 * 7. 理解不同数学函数的定义域和值域
 * 8. 掌握数值方法在实际问题中的应用
 * 
 * 主要内容：
 * - 基础数学函数：三角函数、反三角函数
 * - 指数与对数函数：自然对数、常用对数、指数函数
 * - 幂函数与根函数：平方根、立方根、任意次幂
 * - 双曲函数：双曲正弦、余弦、正切及其反函数
 * - 特殊值处理：NaN传播、无穷大运算、零值边界
 * - 数学常数：π、e、黄金比例等的精确使用
 * - 数值计算方法：牛顿法、数值积分、函数逼近
 * - 实际应用：物理计算、工程数学、统计分析
 */

fn main() {
    println!("=== Rust 浮点运算：数学函数与特殊值 ===\n");

    // 1. 基础数学函数
    demonstrate_basic_math_functions();
    
    // 2. 三角函数详解
    demonstrate_trigonometric_functions();
    
    // 3. 指数与对数函数
    demonstrate_exponential_logarithmic();
    
    // 4. 幂函数与根函数
    demonstrate_power_root_functions();
    
    // 5. 双曲函数
    demonstrate_hyperbolic_functions();
    
    // 6. 特殊值的产生与处理
    demonstrate_special_values_handling();
    
    // 7. 数学常数的精确使用
    demonstrate_mathematical_constants();
    
    // 8. 数值计算方法
    demonstrate_numerical_methods();
    
    // 9. 实际应用案例
    demonstrate_practical_applications();
    
    println!("\n=== 总结 ===");
    println!("浮点数学函数是科学计算和工程应用的基础，掌握其特性对于：");
    println!("- 物理仿真和数学建模");
    println!("- 信号处理和图像分析");
    println!("- 机器学习和数据科学");
    println!("- 金融建模和风险分析");
    println!("都具有重要意义。正确处理特殊值和边界情况是编写健壮程序的关键。");
}

/**
 * 演示基础数学函数的使用
 * 
 * 重点知识：
 * - Rust标准库中的数学函数
 * - 函数的定义域和值域
 * - 精度和性能考虑
 * - 错误处理和边界情况
 */
fn demonstrate_basic_math_functions() {
    println!("1. 基础数学函数");
    println!("   探索Rust标准库中丰富的数学函数");
    
    // 1.1 绝对值函数
    println!("\n1.1 绝对值函数");
    
    let values = [-5.7, -0.0, 0.0, 3.14, f64::INFINITY, f64::NEG_INFINITY];
    println!("绝对值函数 abs():");
    for val in values {
        println!("  |{}| = {}", val, val.abs());
    }
    
    // 处理特殊值
    let nan_abs = f64::NAN.abs();
    println!("  |NaN| = {} (NaN的绝对值仍是NaN)", nan_abs);
    
    // 1.2 舍入函数族
    println!("\n1.2 舍入函数族");
    
    let test_values = [3.2f64, 3.5, 3.8, -2.3, -2.5, -2.7, 0.0, -0.0];
    
    println!("舍入函数对比:");
    println!("值    | floor | ceil  | round | trunc | fract");
    println!("------|-------|-------|-------|-------|-------");
    
    for val in test_values {
        println!("{:5.1} | {:5.1} | {:5.1} | {:5.1} | {:5.1} | {:5.2}",
                val,
                val.floor(),    // 向下取整（向负无穷方向）
                val.ceil(),     // 向上取整（向正无穷方向）
                val.round(),    // 四舍五入到最近整数
                val.trunc(),    // 截断（向零方向取整）
                val.fract()     // 小数部分
        );
    }
    
    println!("\n舍入函数说明:");
    println!("  floor(): 不超过x的最大整数（向负无穷取整）");
    println!("  ceil():  不小于x的最小整数（向正无穷取整）");
    println!("  round(): 四舍五入到最近整数（.5时取偶数）");
    println!("  trunc(): 截断小数部分（向零取整）");
    println!("  fract(): 返回小数部分（x - trunc(x)）");
    
    // 1.3 最值函数
    println!("\n1.3 最值函数");
    
    let a = 3.14f64;
    let b = 2.71f64;
    
    println!("最值比较函数:");
    println!("  min({}, {}) = {}", a, b, a.min(b));
    println!("  max({}, {}) = {}", a, b, a.max(b));
    
    // 特殊值的最值比较
    let special_comparisons = [
        (f64::NAN, 1.0, "NaN vs 1.0"),
        (f64::INFINITY, 1000.0, "∞ vs 1000.0"),
        (f64::NEG_INFINITY, -1000.0, "-∞ vs -1000.0"),
        (0.0, -0.0, "0.0 vs -0.0"),
    ];
    
    println!("特殊值的最值比较:");
    for (x, y, desc) in special_comparisons {
        println!("  {} → min: {}, max: {}", desc, x.min(y), x.max(y));
    }
    
    // 1.4 符号和分类函数
    println!("\n1.4 符号和分类函数");
    
    let classification_values = [
        ("正数", 42.0),
        ("负数", -3.14),
        ("正零", 0.0),
        ("负零", -0.0),
        ("正无穷", f64::INFINITY),
        ("负无穷", f64::NEG_INFINITY),
        ("NaN", f64::NAN),
        ("极小正数", f64::MIN_POSITIVE),
    ];
    
    println!("数值分类和符号判断:");
    println!("值         | 符号+ | 符号- | 有限 | 无限 | NaN  | 正常 | 次正常");
    println!("-----------|-------|-------|------|------|------|------|--------");
    
    for (name, val) in classification_values {
        println!("{:10} | {:5} | {:5} | {:4} | {:4} | {:4} | {:4} | {:6}",
                name,
                val.is_sign_positive(),
                val.is_sign_negative(),
                val.is_finite(),
                val.is_infinite(),
                val.is_nan(),
                val.is_normal(),
                val.is_subnormal(),
        );
    }
    
    // 1.5 符号操作函数
    println!("\n1.5 符号操作函数");
    
    let base = 5.5f64;
    println!("符号操作示例 (基准值: {}):", base);
    println!("  copysign({}, 1.0) = {} (复制正号)", base, base.copysign(1.0));
    println!("  copysign({}, -1.0) = {} (复制负号)", base, base.copysign(-1.0));
    println!("  copysign({}, -0.0) = {} (复制负零的符号)", base, base.copysign(-0.0));
    
    // signum 函数：符号函数
    let signum_values = [5.0, -3.0, 0.0, -0.0, f64::INFINITY, f64::NEG_INFINITY, f64::NAN];
    println!("符号函数 signum():");
    for val in signum_values {
        println!("  signum({}) = {}", val, val.signum());
    }
    
    println!("\n{}", "=".repeat(80));
}

/**
 * 演示三角函数的详细使用
 * 
 * 重点知识：
 * - 弧度制与角度制的转换
 * - 三角函数的周期性和对称性
 * - 反三角函数的值域限制
 * - 三角恒等式的数值验证
 */
fn demonstrate_trigonometric_functions() {
    println!("\n2. 三角函数详解");
    println!("   深入理解三角函数的数学特性和数值计算");
    
    // 2.1 角度制与弧度制
    println!("\n2.1 角度制与弧度制转换");
    
    // 常用角度的弧度值
    let angles_degrees = [0.0f64, 30.0, 45.0, 60.0, 90.0, 120.0, 180.0, 270.0, 360.0];
    
    println!("角度制到弧度制转换:");
    println!("度数 | 弧度       | sin        | cos        | tan");
    println!("-----|------------|------------|------------|------------");
    
    for deg in angles_degrees {
        let rad = deg.to_radians();
        let sin_val = rad.sin();
        let cos_val = rad.cos();
        let tan_val = rad.tan();
        
        println!("{:4.0}° | {:10.6} | {:10.6} | {:10.6} | {:10.6}",
                deg, rad, sin_val, cos_val, tan_val);
    }
    
    // 2.2 三角函数的特殊值
    println!("\n2.2 三角函数的特殊值和精度问题");
    
    use std::f64::consts::PI;
    
    let special_angles = [
        ("0", 0.0),
        ("π/6", PI / 6.0),
        ("π/4", PI / 4.0),
        ("π/3", PI / 3.0),
        ("π/2", PI / 2.0),
        ("π", PI),
        ("3π/2", 3.0 * PI / 2.0),
        ("2π", 2.0 * PI),
    ];
    
    println!("特殊角的三角函数值:");
    for (name, angle) in special_angles {
        println!("{}:", name);
        println!("  sin({}) = {:.15}", name, angle.sin());
        println!("  cos({}) = {:.15}", name, angle.cos());
        println!("  tan({}) = {:.15}", name, angle.tan());
        
        // 检查理论零值
        let sin_val = angle.sin();
        let cos_val = angle.cos();
        if sin_val.abs() < 1e-15 {
            println!("    → sin({}) 理论上应该为 0", name);
        }
        if cos_val.abs() < 1e-15 {
            println!("    → cos({}) 理论上应该为 0", name);
        }
        println!();
    }
    
    // 2.3 三角恒等式的数值验证
    println!("\n2.3 三角恒等式的数值验证");
    
    let test_angles = [PI/6.0, PI/4.0, PI/3.0, PI/2.0, 2.5, -1.2];
    
    println!("验证基本三角恒等式:");
    for angle in test_angles {
        let sin_val = angle.sin();
        let cos_val = angle.cos();
        let tan_val = angle.tan();
        
        // sin²x + cos²x = 1
        let identity1 = sin_val * sin_val + cos_val * cos_val;
        
        // tan(x) = sin(x) / cos(x) (当cos(x) ≠ 0时)
        let identity2 = if cos_val.abs() > 1e-10 {
            sin_val / cos_val
        } else {
            f64::NAN
        };
        
        println!("角度 {:.3} 弧度:", angle);
        println!("  sin²x + cos²x = {:.15} (应该 = 1.0)", identity1);
        println!("  sin/cos = {:.10}, tan = {:.10}", identity2, tan_val);
        println!("  误差: sin²+cos²-1 = {:.2e}", identity1 - 1.0);
        println!();
    }
    
    // 2.4 反三角函数
    println!("\n2.4 反三角函数及其定义域值域");
    
    let test_values = [-1.0f64, -0.866, -0.5, 0.0, 0.5, 0.866, 1.0];
    
    println!("反三角函数计算:");
    println!("值      | arcsin(弧度) | arcsin(度) | arccos(弧度) | arccos(度) | arctan(弧度) | arctan(度)");
    println!("--------|-------------|-----------|-------------|-----------|-------------|----------");
    
    for val in test_values {
        let asin_rad = val.asin();
        let asin_deg = asin_rad.to_degrees();
        let acos_rad = val.acos();
        let acos_deg = acos_rad.to_degrees();
        let atan_rad = val.atan();
        let atan_deg = atan_rad.to_degrees();
        
        println!("{:7.3} | {:11.6} | {:9.3}° | {:11.6} | {:9.3}° | {:11.6} | {:9.3}°",
                val, asin_rad, asin_deg, acos_rad, acos_deg, atan_rad, atan_deg);
    }
    
    // 反三角函数的定义域问题
    println!("\n反三角函数定义域测试:");
    let out_of_domain = [1.5f64, -1.5, 2.0, -2.0];
    
    for val in out_of_domain {
        println!("  arcsin({}) = {} (定义域外)", val, val.asin());
        println!("  arccos({}) = {} (定义域外)", val, val.acos());
    }
    
    // 2.5 atan2函数：四象限反正切
    println!("\n2.5 atan2函数：四象限反正切");
    
    let coordinates = [
        (1.0f64, 1.0, "第一象限"),
        (-1.0, 1.0, "第二象限"),
        (-1.0, -1.0, "第三象限"),
        (1.0, -1.0, "第四象限"),
        (0.0, 1.0, "正y轴"),
        (0.0, -1.0, "负y轴"),
        (1.0, 0.0, "正x轴"),
        (-1.0, 0.0, "负x轴"),
    ];
    
    println!("atan2(y, x) 计算角度 (从正x轴逆时针测量):");
    for (y, x, desc) in coordinates {
        let angle_rad = y.atan2(x);
        let angle_deg = angle_rad.to_degrees();
        println!("  atan2({}, {}) = {:.6} 弧度 = {:.2}° ({})", 
                 y, x, angle_rad, angle_deg, desc);
    }
    
    // 2.6 三角函数的周期性
    println!("\n2.6 三角函数的周期性验证");
    
    let base_angle = PI / 3.0;  // 60度
    println!("验证三角函数周期性 (基准角: π/3):");
    
    for k in 0..3 {
        let angle_with_period = base_angle + k as f64 * 2.0 * PI;
        println!("  角度: π/3 + {}×2π = {:.6}", k, angle_with_period);
        println!("    sin = {:.10}", angle_with_period.sin());
        println!("    cos = {:.10}", angle_with_period.cos());
        println!("    与基准值的差异: sin差={:.2e}, cos差={:.2e}", 
                 (angle_with_period.sin() - base_angle.sin()).abs(),
                 (angle_with_period.cos() - base_angle.cos()).abs());
        println!();
    }
    
    println!("\n{}", "=".repeat(80));
}

/**
 * 演示指数与对数函数
 * 
 * 重点知识：
 * - 自然指数函数和自然对数
 * - 不同底的对数函数
 * - 指数对数的特殊性质
 * - 数值稳定性问题
 */
fn demonstrate_exponential_logarithmic() {
    println!("\n3. 指数与对数函数");
    println!("   掌握指数对数函数的数学性质和数值计算");
    
    // 3.1 自然指数函数
    println!("\n3.1 自然指数函数 exp(x)");
    
    use std::f64::consts::{E, LN_2, LN_10};
    
    let exp_values = [-2.0f64, -1.0, -0.5, 0.0, 0.5, 1.0, 2.0, 5.0, 10.0];
    
    println!("自然指数函数 e^x:");
    println!("x     | e^x        | 近似值");
    println!("------|------------|------------");
    
    for x in exp_values {
        let exp_x = x.exp();
        let approximation = match x {
            0.0 => "1 (定义)".to_string(),
            1.0 => format!("e ≈ {:.6}", E),
            -1.0 => format!("1/e ≈ {:.6}", 1.0/E),
            2.0 => format!("e² ≈ {:.6}", E*E),
            _ => "数值计算".to_string(),
        };
        println!("{:5.1} | {:10.6} | {}", x, exp_x, approximation);
    }
    
    // 指数函数的特殊性质
    println!("\n指数函数性质验证:");
    let x1: f64 = 2.0;
    let x2: f64 = 3.0;
    println!("  e^({} + {}) = e^{} = {:.6}", x1, x2, x1+x2, (x1+x2).exp());
    println!("  e^{} × e^{} = {:.6} × {:.6} = {:.6}", x1, x2, x1.exp(), x2.exp(), x1.exp() * x2.exp());
    println!("  验证: e^(a+b) = e^a × e^b, 误差: {:.2e}", 
             ((x1+x2).exp() - x1.exp() * x2.exp()).abs());
    
    // 3.2 自然对数函数
    println!("\n3.2 自然对数函数 ln(x)");
    
    let ln_values = [0.01, 0.1, 0.5, 1.0, E, 2.0, 10.0, 100.0, 1000.0];
    
    println!("自然对数函数 ln(x):");
    println!("x        | ln(x)      | 特殊意义");
    println!("---------|------------|------------------------");
    
    for x in ln_values {
        let ln_x = x.ln();
        let meaning = match x {
            val if (val - 1.0).abs() < 1e-10 => "ln(1) = 0".to_string(),
            val if (val - E).abs() < 1e-10 => "ln(e) = 1".to_string(),
            val if (val - 10.0).abs() < 1e-10 => format!("ln(10) ≈ {:.6}", LN_10),
            val if (val - 2.0).abs() < 1e-10 => format!("ln(2) ≈ {:.6}", LN_2),
            _ => "数值计算".to_string(),
        };
        println!("{:8.2} | {:10.6} | {}", x, ln_x, meaning);
    }
    
    // 对数函数的定义域问题
    println!("\n对数函数定义域测试:");
    let invalid_ln_inputs = [0.0, -1.0, -10.0, f64::NEG_INFINITY];
    for x in invalid_ln_inputs {
        println!("  ln({}) = {} (定义域外或边界)", x, x.ln());
    }
    
    // 3.3 其他底的对数
    println!("\n3.3 其他底的对数函数");
    
    let log_test_values = [1.0f64, 2.0, 8.0, 10.0, 100.0, 1024.0];
    
    println!("不同底的对数:");
    println!("x       | log₂(x)    | log₁₀(x)   | ln(x)");
    println!("--------|------------|------------|------------");
    
    for x in log_test_values {
        let log2_x = x.log2();
        let log10_x = x.log10();
        let ln_x = x.ln();
        println!("{:7.0} | {:10.6} | {:10.6} | {:10.6}", x, log2_x, log10_x, ln_x);
    }
    
    // 验证换底公式
    println!("\n换底公式验证: log_a(x) = ln(x) / ln(a)");
    let base: f64 = 3.0;
    let x: f64 = 27.0;
    let log_base_x_direct = x.ln() / base.ln();  // 换底公式
    let expected = 3.0;  // 因为 3³ = 27
    
    println!("  log₃(27) 通过换底公式: {:.10}", log_base_x_direct);
    println!("  理论值 (3³ = 27): {:.10}", expected);
    println!("  误差: {:.2e}", (log_base_x_direct - expected).abs());
    
    // 3.4 指数对数的互逆性
    println!("\n3.4 指数对数的互逆性验证");
    
    let inverse_test_values = [0.5f64, 1.0, 2.5, 10.0, -2.0];
    
    println!("指数对数互逆性测试:");
    for x in inverse_test_values {
        if x > 0.0 {
            let exp_ln_x = x.ln().exp();
            let ln_exp_x = x.exp().ln();
            println!("  x = {:.1}:", x);
            println!("    e^(ln(x)) = e^(ln({})) = {:.10} (应该 = {})", x, exp_ln_x, x);
            println!("    ln(e^x) = ln(e^{}) = {:.10} (应该 = {})", x, ln_exp_x, x);
            println!("    误差: e^ln-x = {:.2e}, ln_e^-x = {:.2e}", 
                     (exp_ln_x - x).abs(), (ln_exp_x - x).abs());
        } else {
            let ln_exp_x = x.exp().ln();
            println!("  x = {:.1}: ln(e^x) = {:.10} (应该 = {})", x, ln_exp_x, x);
        }
        println!();
    }
    
    // 3.5 数值稳定性问题
    println!("\n3.5 数值稳定性问题");
    
    // exp(x) - 1 在 x 接近 0 时的精度问题
    println!("exp(x) - 1 在小x值时的精度:");
    let small_values = [1e-1f64, 1e-5, 1e-10, 1e-15];
    
    for x in small_values {
        let exp_minus_1_direct = x.exp() - 1.0;
        let exp_minus_1_stable = x.exp_m1();  // 数值稳定的实现
        let taylor_approximation = x + x*x/2.0;  // 泰勒展开近似
        
        println!("  x = {:.0e}:", x);
        println!("    exp(x) - 1 = {:.15}", exp_minus_1_direct);
        println!("    expm1(x) = {:.15} (稳定算法)", exp_minus_1_stable);
        println!("    泰勒近似 = {:.15} (x + x²/2)", taylor_approximation);
        println!("    相对误差: {:.2e}", 
                 ((exp_minus_1_direct - exp_minus_1_stable) / exp_minus_1_stable).abs());
        println!();
    }
    
    // ln(1 + x) 在 x 接近 0 时的精度问题
    println!("ln(1 + x) 在小x值时的精度:");
    for x in small_values {
        let ln_1_plus_x_direct = (1.0 + x).ln();
        let ln_1_plus_x_stable = x.ln_1p();  // 数值稳定的实现
        let taylor_approximation = x - x*x/2.0;  // 泰勒展开近似
        
        println!("  x = {:.0e}:", x);
        println!("    ln(1 + x) = {:.15}", ln_1_plus_x_direct);
        println!("    ln1p(x) = {:.15} (稳定算法)", ln_1_plus_x_stable);
        println!("    泰勒近似 = {:.15} (x - x²/2)", taylor_approximation);
        println!("    相对误差: {:.2e}", 
                 ((ln_1_plus_x_direct - ln_1_plus_x_stable) / ln_1_plus_x_stable).abs());
        println!();
    }
    
    println!("\n{}", "=".repeat(80));
}

/**
 * 演示幂函数与根函数
 * 
 * 重点知识：
 * - 整数幂和有理数幂
 * - 平方根、立方根、n次根
 * - 负数的根函数问题
 * - 精度和性能考虑
 */
fn demonstrate_power_root_functions() {
    println!("\n4. 幂函数与根函数");
    println!("   理解各种幂运算和根运算的数学特性");
    
    // 4.1 幂函数 powf, powi, powf
    println!("\n4.1 幂函数的不同实现");
    
    let base: f64 = 2.0;
    
    // powi: 整数幂（性能最优）
    println!("整数幂函数 powi(整数指数):");
    for i in -3..=5 {
        println!("  {}^{} = {:.6}", base, i, base.powi(i));
    }
    
    // powf: 浮点数幂
    println!("\n浮点数幂函数 powf(浮点指数):");
    let float_exponents = [0.5, 1.5, 2.5, -0.5, 0.0, 1.0];
    for exp in float_exponents {
        println!("  {}^{} = {:.6}", base, exp, base.powf(exp));
    }
    
    // 特殊底数的幂运算
    println!("\n特殊底数的幂运算:");
    let special_bases = [0.0, -1.0, -2.0, f64::INFINITY, f64::NEG_INFINITY];
    let exponent = 2.0;
    
    for base in special_bases {
        println!("  {}^{} = {}", base, exponent, base.powf(exponent));
    }
    
    // 4.2 根函数
    println!("\n4.2 根函数");
    
    let test_values = [0.0f64, 1.0, 4.0, 9.0, 16.0, 25.0, 100.0, 2.0, 8.0, 27.0];
    
    println!("根函数计算:");
    println!("x     | √x        | ∛x        | x^(1/4)   | x^(1/5)");
    println!("------|-----------|-----------|-----------|----------");
    
    for x in test_values {
        let sqrt_x = x.sqrt();                 // 平方根
        let cbrt_x = x.cbrt();                 // 立方根
        let fourth_root = x.powf(0.25);        // 四次根
        let fifth_root = x.powf(0.2);          // 五次根
        
        println!("{:5.1} | {:9.6} | {:9.6} | {:9.6} | {:9.6}",
                x, sqrt_x, cbrt_x, fourth_root, fifth_root);
    }
    
    // 4.3 负数的根运算问题
    println!("\n4.3 负数的根运算");
    
    let negative_values = [-1.0f64, -4.0, -8.0, -27.0];
    
    println!("负数的根运算:");
    for x in negative_values {
        let sqrt_x = x.sqrt();    // 会产生 NaN
        let cbrt_x = x.cbrt();    // 立方根可以计算负数
        let pow_half = x.powf(0.5);  // 等价于sqrt，也是NaN
        let pow_third = x.powf(1.0/3.0);  // 可能不等于cbrt
        
        println!("  x = {}:", x);
        println!("    √x = {} (平方根)", sqrt_x);
        println!("    ∛x = {} (立方根)", cbrt_x);
        println!("    x^0.5 = {} (等价于√x)", pow_half);
        println!("    x^(1/3) = {} (可能与∛x不同)", pow_third);
        println!();
    }
    
    // 4.4 根函数与幂函数的关系验证
    println!("\n4.4 根函数与幂函数的关系验证");
    
    let test_bases = [4.0f64, 9.0, 16.0, 25.0];
    
    for base in test_bases {
        let sqrt_val = base.sqrt();
        let pow_half = base.powf(0.5);
        let cbrt_val = base.cbrt();
        let pow_third = base.powf(1.0/3.0);
        
        println!("  x = {}:", base);
        println!("    √x = {:.10}, x^0.5 = {:.10}, 差异: {:.2e}", 
                sqrt_val, pow_half, (sqrt_val - pow_half).abs());
        println!("    ∛x = {:.10}, x^(1/3) = {:.10}, 差异: {:.2e}", 
                cbrt_val, pow_third, (cbrt_val - pow_third).abs());
        
        // 验证逆运算
        let sqrt_squared = sqrt_val * sqrt_val;
        let cbrt_cubed = cbrt_val * cbrt_val * cbrt_val;
        
        println!("    (√x)² = {:.10}, 还原误差: {:.2e}", 
                sqrt_squared, (sqrt_squared - base).abs());
        println!("    (∛x)³ = {:.10}, 还原误差: {:.2e}", 
                cbrt_cubed, (cbrt_cubed - base).abs());
        println!();
    }
    
    // 4.5 幂函数的特殊性质
    println!("\n4.5 幂函数的特殊性质验证");
    
    let a: f64 = 3.0;
    let b: f64 = 4.0;
    let m: f64 = 2.0;
    let n: f64 = 1.5;
    
    println!("验证幂函数运算法则:");
    
    // a^m × a^n = a^(m+n)
    let left1 = a.powf(m) * a.powf(n);
    let right1 = a.powf(m + n);
    println!("  a^m × a^n = a^(m+n)");
    println!("  {}^{} × {}^{} = {} × {} = {:.6}", a, m, a, n, a.powf(m), a.powf(n), left1);
    println!("  {}^({} + {}) = {}^{} = {:.6}", a, m, n, a, m+n, right1);
    println!("  误差: {:.2e}", (left1 - right1).abs());
    println!();
    
    // (a^m)^n = a^(m×n)
    let left2 = a.powf(m).powf(n);
    let right2 = a.powf(m * n);
    println!("  (a^m)^n = a^(m×n)");
    println!("  ({}^{})^{} = {:.6}^{} = {:.6}", a, m, n, a.powf(m), n, left2);
    println!("  {}^({} × {}) = {}^{} = {:.6}", a, m, n, a, m*n, right2);
    println!("  误差: {:.2e}", (left2 - right2).abs());
    println!();
    
    // (ab)^m = a^m × b^m
    let left3 = (a * b).powf(m);
    let right3 = a.powf(m) * b.powf(m);
    println!("  (ab)^m = a^m × b^m");
    println!("  ({} × {})^{} = {}^{} = {:.6}", a, b, m, a*b, m, left3);
    println!("  {}^{} × {}^{} = {} × {} = {:.6}", a, m, b, m, a.powf(m), b.powf(m), right3);
    println!("  误差: {:.2e}", (left3 - right3).abs());
    
    println!("\n{}", "=".repeat(80));
}

/**
 * 演示双曲函数
 * 
 * 重点知识：
 * - 双曲正弦、余弦、正切函数
 * - 反双曲函数
 * - 双曲函数的性质和应用
 * - 与三角函数的关系
 */
fn demonstrate_hyperbolic_functions() {
    println!("\n5. 双曲函数");
    println!("   探索双曲函数的数学特性和实际应用");
    
    // 5.1 双曲函数的基本概念
    println!("\n5.1 双曲函数的定义和基本性质");
    
    
    let test_values = [-2.0f64, -1.0, -0.5, 0.0, 0.5, 1.0, 2.0];
    
    println!("双曲函数值:");
    println!("x     | sinh(x)    | cosh(x)    | tanh(x)    | 定义验证");
    println!("------|------------|------------|------------|------------");
    
    for x in test_values {
        let sinh_x = x.sinh();
        let cosh_x = x.cosh();
        let tanh_x = x.tanh();
        
        // 根据定义验证
        let sinh_def = (x.exp() - (-x).exp()) / 2.0;
        let cosh_def = (x.exp() + (-x).exp()) / 2.0;
        let tanh_def = sinh_x / cosh_x;
        
        println!("{:5.1} | {:10.6} | {:10.6} | {:10.6} | sinh定义: {:.6}",
                x, sinh_x, cosh_x, tanh_x, sinh_def);
        println!("      |            |            |            | cosh定义: {:.6}", cosh_def);
        println!("      |            |            |            | tanh=sinh/cosh: {:.6}", tanh_def);
        println!();
    }
    
    // 5.2 双曲函数的重要恒等式
    println!("\n5.2 双曲函数恒等式验证");
    
    let identity_test_values = [0.5f64, 1.0, 1.5, 2.0];
    
    for x in identity_test_values {
        let sinh_x = x.sinh();
        let cosh_x = x.cosh();
        let tanh_x = x.tanh();
        
        // cosh²x - sinh²x = 1 (双曲函数基本恒等式)
        let identity1 = cosh_x * cosh_x - sinh_x * sinh_x;
        
        // tanh²x + sech²x = 1, 即 tanh²x + 1/cosh²x = 1
        let identity2 = tanh_x * tanh_x + 1.0 / (cosh_x * cosh_x);
        
        // sinh(2x) = 2sinh(x)cosh(x)
        let sinh_2x = (2.0 * x).sinh();
        let sinh_double_angle = 2.0 * sinh_x * cosh_x;
        
        // cosh(2x) = cosh²x + sinh²x
        let cosh_2x = (2.0 * x).cosh();
        let cosh_double_angle = cosh_x * cosh_x + sinh_x * sinh_x;
        
        println!("x = {:.1} 的双曲恒等式验证:", x);
        println!("  cosh²x - sinh²x = {:.15} (应该 = 1.0)", identity1);
        println!("  tanh²x + sech²x = {:.15} (应该 = 1.0)", identity2);
        println!("  sinh(2x) = {:.10}, 2sinh(x)cosh(x) = {:.10}", sinh_2x, sinh_double_angle);
        println!("  cosh(2x) = {:.10}, cosh²x + sinh²x = {:.10}", cosh_2x, cosh_double_angle);
        println!("  误差: cosh²-sinh²-1 = {:.2e}", (identity1 - 1.0).abs());
        println!("        sinh倍角误差 = {:.2e}", (sinh_2x - sinh_double_angle).abs());
        println!("        cosh倍角误差 = {:.2e}", (cosh_2x - cosh_double_angle).abs());
        println!();
    }
    
    // 5.3 反双曲函数
    println!("\n5.3 反双曲函数");
    
    let inverse_test_values = [-2.0f64, -1.0, -0.5, 0.0, 0.5, 1.0, 2.0];
    
    println!("反双曲函数:");
    println!("x     | asinh(x)   | acosh(x)   | atanh(x)   | 定义域说明");
    println!("------|------------|------------|------------|-------------");
    
    for x in inverse_test_values {
        let asinh_x = x.asinh();  // 定义域: 全实数
        let acosh_x = x.acosh();  // 定义域: x ≥ 1
        let atanh_x = x.atanh();  // 定义域: |x| < 1
        
        let domain_note = if x >= 1.0 {
            "全部有定义"
        } else if x > -1.0 && x < 1.0 {
            "asinh, atanh有定义"
        } else {
            "仅asinh有定义"
        };
        
        println!("{:5.1} | {:10.6} | {:10.6} | {:10.6} | {}",
                x, asinh_x, acosh_x, atanh_x, domain_note);
    }
    
    // 反双曲函数的边界值测试
    println!("\n反双曲函数边界值测试:");
    let boundary_values = [0.999f64, 1.0, 1.001, -0.999, -1.0, -1.001];
    
    for x in boundary_values {
        println!("  x = {:.3}:", x);
        println!("    atanh(x) = {} (定义域: |x| < 1)", x.atanh());
        println!("    acosh(x) = {} (定义域: x ≥ 1)", x.acosh());
    }
    
    // 5.4 双曲函数与反双曲函数的互逆性
    println!("\n5.4 双曲函数与反双曲函数的互逆性验证");
    
    let inverse_verify_values = [0.5f64, 1.0, 1.5, 2.0];
    
    for x in inverse_verify_values {
        println!("x = {:.1}:", x);
        
        // sinh和asinh的互逆性
        let sinh_asinh = x.sinh().asinh();
        let asinh_sinh = x.asinh().sinh();
        println!("  asinh(sinh(x)) = {:.10} (应该 = {})", sinh_asinh, x);
        println!("  sinh(asinh(x)) = {:.10} (应该 = {})", asinh_sinh, x);
        
        // cosh和acosh的互逆性（仅对x≥1）
        if x >= 1.0 {
            let cosh_acosh = x.cosh().acosh();
            let acosh_cosh = x.acosh().cosh();
            println!("  acosh(cosh(x)) = {:.10} (应该 = {})", cosh_acosh, x);
            println!("  cosh(acosh(x)) = {:.10} (应该 = {})", acosh_cosh, x);
        }
        
        // tanh和atanh的互逆性（仅对|x|<1）
        if x.abs() < 1.0 {
            let tanh_atanh = x.tanh().atanh();
            let atanh_tanh = x.atanh().tanh();
            println!("  atanh(tanh(x)) = {:.10} (应该 = {})", tanh_atanh, x);
            println!("  tanh(atanh(x)) = {:.10} (应该 = {})", atanh_tanh, x);
        }
        println!();
    }
    
    // 5.5 双曲函数的实际应用
    println!("\n5.5 双曲函数的实际应用");
    
    // 悬链线 (Catenary) - 绳子在重力作用下的形状
    println!("悬链线方程: y = a·cosh(x/a)");
    let a = 10.0;  // 悬链线参数
    let x_positions = [-20.0f64, -10.0, 0.0, 10.0, 20.0];
    
    println!("悬链线上的点 (a = {}):", a);
    for x in x_positions {
        let y = a * (x / a).cosh();
        println!("  x = {:5.1}, y = {:8.3}", x, y);
    }
    
    // 双曲函数在相对论中的应用
    println!("\n双曲函数在物理学中的应用:");
    let c = 299792458.0; // 光速 (m/s)
    let velocities = [0.1 * c, 0.5 * c, 0.9 * c, 0.99 * c];
    
    println!("相对论洛伦兹因子 γ = 1/√(1-v²/c²):");
    for v in velocities {
        let beta: f64 = v / c;  // v/c
        let gamma = 1.0 / (1.0 - beta * beta).sqrt();
        let rapidity = beta.atanh();  // 快度参数
        
        println!("  v = {:.2}c, β = {:.2}, γ = {:.4}, φ = {:.4}", 
                v/c, beta, gamma, rapidity);
    }
    
    println!("\n{}", "=".repeat(80));
}

/**
 * 演示特殊值的产生与处理
 * 
 * 重点知识：
 * - NaN的传播和检测
 * - 无穷大的运算规则
 * - 零值的特殊情况
 * - 特殊值的安全处理策略
 */
fn demonstrate_special_values_handling() {
    println!("\n6. 特殊值的产生与处理");
    println!("   深入理解IEEE 754特殊值的行为和处理方法");
    
    // 6.1 NaN的产生原因
    println!("\n6.1 NaN (Not a Number) 的产生原因");
    
    let nan_operations = [
        ("0.0 / 0.0", 0.0 / 0.0),
        ("∞ - ∞", f64::INFINITY - f64::INFINITY),
        ("∞ / ∞", f64::INFINITY / f64::INFINITY),
        ("∞ × 0", f64::INFINITY * 0.0),
        ("0 × ∞", 0.0 * f64::INFINITY),
        ("(-1.0).sqrt()", (-1.0_f64).sqrt()),
        ("(-1.0).ln()", (-1.0_f64).ln()),
        ("f64::NAN + 1.0", f64::NAN + 1.0),
        ("1.0.powf(f64::NAN)", 1.0_f64.powf(f64::NAN)),
        ("f64::NAN.sin()", f64::NAN.sin()),
    ];
    
    println!("产生NaN的运算:");
    for (desc, result) in nan_operations {
        println!("  {} = {}", desc, result);
        println!("    is_nan(): {}", result.is_nan());
    }
    
    // 6.2 NaN的特殊性质
    println!("\n6.2 NaN的特殊性质");
    
    let nan1 = f64::NAN;
    let nan2 = 0.0 / 0.0;
    let normal = 42.0;
    
    println!("NaN比较的特殊性:");
    println!("  NaN == NaN: {} (总是false!)", nan1 == nan2);
    println!("  NaN != NaN: {} (总是true!)", nan1 != nan2);
    println!("  NaN < 42.0: {} (总是false)", nan1 < normal);
    println!("  NaN > 42.0: {} (总是false)", nan1 > normal);
    println!("  NaN == 42.0: {} (总是false)", nan1 == normal);
    println!("  NaN != 42.0: {} (总是true)", nan1 != normal);
    
    // NaN在排序中的行为
    let mut values_with_nan = [3.0, f64::NAN, 1.0, 2.0, f64::NAN];
    println!("\n排序前: {:?}", values_with_nan);
    values_with_nan.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
    println!("排序后: {:?}", values_with_nan);
    
    // 6.3 无穷大的运算规则
    println!("\n6.3 无穷大的运算规则");
    
    let pos_inf = f64::INFINITY;
    let neg_inf = f64::NEG_INFINITY;
    let finite = 42.0;
    
    println!("无穷大的算术运算:");
    println!("  ∞ + 42 = {}", pos_inf + finite);
    println!("  ∞ - 42 = {}", pos_inf - finite);
    println!("  ∞ × 42 = {}", pos_inf * finite);
    println!("  ∞ ÷ 42 = {}", pos_inf / finite);
    println!("  42 ÷ ∞ = {}", finite / pos_inf);
    println!("  ∞ + ∞ = {}", pos_inf + pos_inf);
    println!("  ∞ - (-∞) = {}", pos_inf - neg_inf);
    println!("  ∞ × ∞ = {}", pos_inf * pos_inf);
    println!("  ∞ × (-∞) = {}", pos_inf * neg_inf);
    
    println!("\n产生NaN的无穷大运算:");
    println!("  ∞ - ∞ = {}", pos_inf - pos_inf);
    println!("  ∞ ÷ ∞ = {}", pos_inf / pos_inf);
    println!("  ∞ × 0 = {}", pos_inf * 0.0);
    println!("  0 × ∞ = {}", 0.0 * pos_inf);
    
    // 6.4 无穷大的数学函数
    println!("\n6.4 无穷大在数学函数中的行为");
    
    let inf_function_tests = [
        ("sin(∞)", pos_inf.sin()),
        ("cos(∞)", pos_inf.cos()),
        ("tan(∞)", pos_inf.tan()),
        ("ln(∞)", pos_inf.ln()),
        ("exp(∞)", pos_inf.exp()),
        ("exp(-∞)", neg_inf.exp()),
        ("sqrt(∞)", pos_inf.sqrt()),
        ("∞^2", pos_inf.powf(2.0)),
        ("∞^0", pos_inf.powf(0.0)),
        ("1^∞", 1.0_f64.powf(pos_inf)),
    ];
    
    for (desc, result) in inf_function_tests {
        println!("  {} = {}", desc, result);
    }
    
    // 6.5 零值的特殊情况
    println!("\n6.5 零值的特殊情况");
    
    let pos_zero: f64 = 0.0;
    let neg_zero: f64 = -0.0;
    
    println!("正零和负零:");
    println!("  +0.0 == -0.0: {}", pos_zero == neg_zero);
    println!("  1.0 / +0.0 = {}", 1.0 / pos_zero);
    println!("  1.0 / -0.0 = {}", 1.0 / neg_zero);
    println!("  (+0.0).copysign(1.0) = {}", pos_zero.copysign(1.0));
    println!("  (-0.0).copysign(1.0) = {}", neg_zero.copysign(1.0));
    println!("  (+0.0).is_sign_positive() = {}", pos_zero.is_sign_positive());
    println!("  (-0.0).is_sign_positive() = {}", neg_zero.is_sign_positive());
    
    // 零在数学函数中的边界行为
    println!("\n零值的数学函数边界行为:");
    let zero_function_tests = [
        ("ln(0.0)", 0.0_f64.ln()),
        ("ln(-0.0)", (-0.0_f64).ln()),
        ("sqrt(0.0)", 0.0_f64.sqrt()),
        ("sqrt(-0.0)", (-0.0_f64).sqrt()),
        ("0.0.powf(0.0)", 0.0_f64.powf(0.0)),
        ("atan2(0.0, 0.0)", 0.0_f64.atan2(0.0)),
        ("atan2(-0.0, 0.0)", (-0.0_f64).atan2(0.0)),
        ("atan2(0.0, -0.0)", 0.0_f64.atan2(-0.0)),
    ];
    
    for (desc, result) in zero_function_tests {
        println!("  {} = {}", desc, result);
    }
    
    // 6.6 特殊值的安全处理策略
    println!("\n6.6 特殊值的安全处理策略");
    
    // 安全除法函数
    fn safe_divide(a: f64, b: f64) -> Option<f64> {
        if b == 0.0 {
            None
        } else if a.is_finite() && b.is_finite() {
            Some(a / b)
        } else if a.is_nan() || b.is_nan() {
            None
        } else {
            Some(a / b)  // 处理无穷大情况
        }
    }
    
    // 安全平方根函数
    fn safe_sqrt(x: f64) -> Result<f64, &'static str> {
        if x.is_nan() {
            Err("输入为NaN")
        } else if x < 0.0 {
            Err("负数没有实数平方根")
        } else {
            Ok(x.sqrt())
        }
    }
    
    // 安全对数函数
    fn safe_ln(x: f64) -> Result<f64, &'static str> {
        if x.is_nan() {
            Err("输入为NaN")
        } else if x <= 0.0 {
            Err("对数的参数必须为正数")
        } else if x.is_infinite() {
            Ok(f64::INFINITY)
        } else {
            Ok(x.ln())
        }
    }
    
    println!("安全函数示例:");
    let test_cases = [1.0, 0.0, -1.0, f64::INFINITY, f64::NAN];
    
    for x in test_cases {
        println!("  x = {}:", x);
        
        match safe_divide(10.0, x) {
            Some(result) => println!("    safe_divide(10.0, x) = {}", result),
            None => println!("    safe_divide(10.0, x) = 无法计算"),
        }
        
        match safe_sqrt(x) {
            Ok(result) => println!("    safe_sqrt(x) = {}", result),
            Err(msg) => println!("    safe_sqrt(x) = 错误: {}", msg),
        }
        
        match safe_ln(x) {
            Ok(result) => println!("    safe_ln(x) = {}", result),
            Err(msg) => println!("    safe_ln(x) = 错误: {}", msg),
        }
        println!();
    }
    
    // 6.7 特殊值传播
    println!("\n6.7 特殊值传播规则");
    
    println!("NaN传播 (任何涉及NaN的运算结果都是NaN):");
    let nan = f64::NAN;
    let operations_with_nan = [
        ("NaN + 1", nan + 1.0),
        ("NaN × 0", nan * 0.0),
        ("NaN.sin()", nan.sin()),
        ("NaN.max(5)", nan.max(5.0)),
        ("5.0.min(NaN)", 5.0_f64.min(nan)),
    ];
    
    for (desc, result) in operations_with_nan {
        println!("  {} = {}", desc, result);
    }
    
    println!("\n无穷大传播 (大部分情况保持无穷大):");
    let inf = f64::INFINITY;
    let inf_operations = [
        ("∞ + 100", inf + 100.0),
        ("∞ × 2", inf * 2.0),
        ("∞.sin()", inf.sin()),  // 这个会产生NaN
        ("∞.exp()", inf.exp()),
        ("(∞).abs()", inf.abs()),
    ];
    
    for (desc, result) in inf_operations {
        println!("  {} = {}", desc, result);
    }
    
    println!("\n{}", "=".repeat(80));
}

/**
 * 演示数学常数的精确使用
 * 
 * 重点知识：
 * - Rust提供的数学常数
 * - 常数的精度和使用场景
 * - 自定义高精度常数
 * - 常数在计算中的重要性
 */
fn demonstrate_mathematical_constants() {
    println!("\n7. 数学常数的精确使用");
    println!("   掌握重要数学常数的精确值和应用");
    
    // 7.1 Rust提供的标准数学常数
    println!("\n7.1 Rust标准库中的数学常数");
    
    use std::f64::consts::*;
    
    let rust_constants = [
        ("π (PI)", PI, "圆周率"),
        ("τ (TAU)", TAU, "2π，完整圆的弧度"),
        ("e (E)", E, "自然常数，欧拉数"),
        ("√2 (SQRT_2)", SQRT_2, "2的平方根"),
        ("1/√2 (FRAC_1_SQRT_2)", FRAC_1_SQRT_2, "1除以√2"),
        ("ln(2) (LN_2)", LN_2, "2的自然对数"),
        ("ln(10) (LN_10)", LN_10, "10的自然对数"),
        ("log₂(e) (LOG2_E)", LOG2_E, "e的以2为底的对数"),
        ("log₁₀(e) (LOG10_E)", LOG10_E, "e的常用对数"),
        ("1/π (FRAC_1_PI)", FRAC_1_PI, "1除以π"),
        ("2/π (FRAC_2_PI)", FRAC_2_PI, "2除以π"),
        ("2/√π (FRAC_2_SQRT_PI)", FRAC_2_SQRT_PI, "2除以√π"),
        ("π/2 (FRAC_PI_2)", FRAC_PI_2, "π除以2"),
        ("π/3 (FRAC_PI_3)", FRAC_PI_3, "π除以3"),
        ("π/4 (FRAC_PI_4)", FRAC_PI_4, "π除以4"),
        ("π/6 (FRAC_PI_6)", FRAC_PI_6, "π除以6"),
        ("π/8 (FRAC_PI_8)", FRAC_PI_8, "π除以8"),
    ];
    
    println!("标准数学常数 (f64精度):");
    for (name, value, description) in rust_constants {
        println!("  {:20} = {:20.15} - {}", name, value, description);
    }
    
    // 7.2 常数精度比较
    println!("\n7.2 常数精度验证和比较");
    
    // 验证一些数学关系
    println!("数学关系验证:");
    
    // π = 4 * arctan(1)
    let pi_from_arctan = 4.0 * 1.0_f64.atan();
    println!("  π 通过 4×arctan(1): {:.15}", pi_from_arctan);
    println!("  标准库 π:           {:.15}", PI);
    println!("  误差:              {:.2e}", (pi_from_arctan - PI).abs());
    
    // e = exp(1)
    let e_from_exp = 1.0_f64.exp();
    println!("  e 通过 exp(1):      {:.15}", e_from_exp);
    println!("  标准库 e:           {:.15}", E);
    println!("  误差:              {:.2e}", (e_from_exp - E).abs());
    
    // ln(2) = ln(2)
    let ln2_computed = 2.0_f64.ln();
    println!("  ln(2) 通过计算:     {:.15}", ln2_computed);
    println!("  标准库 LN_2:        {:.15}", LN_2);
    println!("  误差:              {:.2e}", (ln2_computed - LN_2).abs());
    
    // √2 = sqrt(2)
    let sqrt2_computed = 2.0_f64.sqrt();
    println!("  √2 通过计算:        {:.15}", sqrt2_computed);
    println!("  标准库 SQRT_2:      {:.15}", SQRT_2);
    println!("  误差:              {:.2e}", (sqrt2_computed - SQRT_2).abs());
    
    // 7.3 其他重要数学常数
    println!("\n7.3 其他重要数学常数");
    
    // 黄金比例 φ = (1 + √5) / 2
    let golden_ratio = (1.0 + 5.0_f64.sqrt()) / 2.0;
    println!("黄金比例 φ = (1 + √5) / 2:");
    println!("  φ = {:.15}", golden_ratio);
    println!("  φ² = {:.15}", golden_ratio * golden_ratio);
    println!("  φ + 1 = {:.15} (应该等于φ²)", golden_ratio + 1.0);
    println!("  误差: {:.2e}", ((golden_ratio * golden_ratio) - (golden_ratio + 1.0)).abs());
    
    // 欧拉-马斯刻若尼常数 γ (近似值)
    let euler_mascheroni = 0.5772156649015329;
    println!("\n欧拉-马斯刻若尼常数 γ:");
    println!("  γ ≈ {:.15} (近似值)", euler_mascheroni);
    
    // 阿培里常数 ζ(3)
    let apery_constant = 1.2020569031595943;
    println!("阿培里常数 ζ(3):");
    println!("  ζ(3) ≈ {:.15} (近似值)", apery_constant);
    
    // 卡塔兰常数
    let catalan_constant = 0.9159655941772190;
    println!("卡塔兰常数 G:");
    println!("  G ≈ {:.15} (近似值)", catalan_constant);
    
    // 7.4 物理常数（使用浮点数表示）
    println!("\n7.4 重要物理常数");
    
    let physical_constants = [
        ("光速", 299_792_458.0, "m/s"),
        ("普朗克常数", 6.62607015e-34, "J⋅s"),
        ("玻尔兹曼常数", 1.380649e-23, "J/K"),
        ("阿伏伽德罗常数", 6.02214076e23, "1/mol"),
        ("万有引力常数", 6.67430e-11, "m³/(kg⋅s²)"),
        ("精细结构常数", 7.2973525693e-3, "无量纲"),
        ("电子质量", 9.1093837015e-31, "kg"),
        ("质子质量", 1.67262192369e-27, "kg"),
    ];
    
    println!("物理常数 (SI单位):");
    for (name, value, unit) in physical_constants {
        println!("  {:12} = {:15.6e} {}", name, value, unit);
    }
    
    // 7.5 常数在计算中的应用
    println!("\n7.5 常数在计算中的应用示例");
    
    // 圆的面积和周长计算
    let radius = 5.0;
    let area = PI * radius * radius;
    let circumference = 2.0 * PI * radius;
    println!("圆形计算 (半径 = {}):", radius);
    println!("  面积 = π × r² = {:.6}", area);
    println!("  周长 = 2π × r = {:.6}", circumference);
    
    // 正态分布的概率密度函数
    fn normal_pdf(x: f64, mean: f64, std_dev: f64) -> f64 {
        let coefficient = 1.0 / (std_dev * (2.0 * PI).sqrt());
        let exponent = -0.5 * ((x - mean) / std_dev).powf(2.0);
        coefficient * exponent.exp()
    }
    
    println!("\n正态分布 N(0,1) 的概率密度:");
    let x_values = [-2.0, -1.0, 0.0, 1.0, 2.0];
    for x in x_values {
        let pdf_value = normal_pdf(x, 0.0, 1.0);
        println!("  f({:.1}) = {:.6}", x, pdf_value);
    }
    
    // 复利计算中的 e
    let principal = 1000.0f64;
    let rate = 0.05f64;  // 5%年利率
    let time = 10.0f64;   // 10年
    
    // 连续复利：A = P * e^(rt)
    let continuous_compound = principal * (rate * time).exp();
    println!("\n复利计算 (本金 ¥{}, 年利率 {}%, {} 年):", principal, rate * 100.0, time);
    println!("  连续复利 A = P × e^(rt) = ¥{:.2}", continuous_compound);
    
    // 年复利：A = P * (1 + r)^t
    let annual_compound = principal * (1.0 + rate).powf(time);
    println!("  年复利   A = P × (1+r)^t = ¥{:.2}", annual_compound);
    println!("  差异: ¥{:.2}", continuous_compound - annual_compound);
    
    // 7.6 高精度常数的需要
    println!("\n7.6 高精度需求场景");
    
    // 计算 π 的级数展开（莱布尼兹级数）
    fn compute_pi_leibniz(terms: usize) -> f64 {
        let mut pi_approx = 0.0;
        for n in 0..terms {
            let term = (-1.0_f64).powi(n as i32) / (2 * n + 1) as f64;
            pi_approx += term;
        }
        pi_approx * 4.0
    }
    
    println!("π 的莱布尼兹级数近似:");
    let term_counts = [10, 100, 1000, 10000, 100000];
    for &n in &term_counts {
        let pi_approx = compute_pi_leibniz(n);
        let error = (pi_approx - PI).abs();
        println!("  {} 项: π ≈ {:.10}, 误差: {:.2e}", n, pi_approx, error);
    }
    
    // 计算 e 的级数展开
    fn compute_e_series(terms: usize) -> f64 {
        let mut e_approx = 0.0;
        let mut factorial = 1.0;
        for n in 0..terms {
            if n > 0 {
                factorial *= n as f64;
            }
            e_approx += 1.0 / factorial;
        }
        e_approx
    }
    
    println!("\ne 的级数近似 (e = Σ(1/n!)):");
    for &n in &[5, 10, 15, 20] {
        let e_approx = compute_e_series(n);
        let error = (e_approx - E).abs();
        println!("  {} 项: e ≈ {:.12}, 误差: {:.2e}", n, e_approx, error);
    }
    
    println!("\n数学常数使用建议:");
    println!("  ✓ 优先使用标准库提供的常数（精度高、性能好）");
    println!("  ✓ 对于未提供的常数，预计算并存储高精度值");
    println!("  ✓ 在高精度需求场景中，考虑使用任意精度库");
    println!("  ✓ 注意常数运算可能引入的累积误差");
    
    println!("\n{}", "=".repeat(80));
}

/**
 * 演示数值计算方法
 * 
 * 重点知识：
 * - 牛顿法求根
 * - 数值积分方法
 * - 函数逼近技术
 * - 迭代算法的收敛性
 */
fn demonstrate_numerical_methods() {
    println!("\n8. 数值计算方法");
    println!("   学习重要的数值算法和计算技巧");
    
    // 8.1 牛顿法求根
    println!("\n8.1 牛顿法求根");
    
    // 求解方程 x² - 2 = 0，即求 √2
    fn newton_sqrt2(x0: f64, tolerance: f64, max_iterations: usize) -> (f64, usize) {
        let mut x = x0;
        let mut iterations = 0;
        
        for i in 0..max_iterations {
            let f = x * x - 2.0;  // f(x) = x² - 2
            let f_prime = 2.0 * x;  // f'(x) = 2x
            
            if f.abs() < tolerance {
                return (x, i);
            }
            
            x = x - f / f_prime;  // 牛顿迭代公式
            iterations = i + 1;
        }
        
        (x, iterations)
    }
    
    println!("牛顿法求 √2 (求解 x² - 2 = 0):");
    let (sqrt2_result, iterations) = newton_sqrt2(1.5, 1e-15, 20);
    let actual_sqrt2 = 2.0_f64.sqrt();
    
    println!("  初始值: 1.5");
    println!("  容差: 1e-15");
    println!("  迭代次数: {}", iterations);
    println!("  计算结果: {:.15}", sqrt2_result);
    println!("  标准库值: {:.15}", actual_sqrt2);
    println!("  误差: {:.2e}", (sqrt2_result - actual_sqrt2).abs());
    
    // 求解三次方程 x³ - 2x - 5 = 0
    fn newton_cubic(x0: f64, tolerance: f64, max_iterations: usize) -> (f64, usize) {
        let mut x = x0;
        
        for i in 0..max_iterations {
            let f = x * x * x - 2.0 * x - 5.0;     // f(x) = x³ - 2x - 5
            let f_prime = 3.0 * x * x - 2.0;       // f'(x) = 3x² - 2
            
            if f.abs() < tolerance {
                return (x, i);
            }
            
            if f_prime.abs() < 1e-12 {  // 避免除零
                break;
            }
            
            x = x - f / f_prime;
        }
        
        (x, max_iterations)
    }
    
    println!("\n牛顿法求解 x³ - 2x - 5 = 0:");
    let (root, iterations) = newton_cubic(2.0, 1e-12, 50);
    let verification = root * root * root - 2.0 * root - 5.0;
    
    println!("  初始值: 2.0");
    println!("  迭代次数: {}", iterations);
    println!("  求得的根: {:.12}", root);
    println!("  验证 f(x): {:.2e}", verification);
    
    // 8.2 数值积分
    println!("\n8.2 数值积分方法");
    
    // 矩形法则（左端点）
    fn rectangle_rule<F>(f: F, a: f64, b: f64, n: usize) -> f64 
    where F: Fn(f64) -> f64 {
        let h = (b - a) / n as f64;
        let mut sum = 0.0;
        
        for i in 0..n {
            let x = a + i as f64 * h;
            sum += f(x);
        }
        
        sum * h
    }
    
    // 梯形法则
    fn trapezoidal_rule<F>(f: F, a: f64, b: f64, n: usize) -> f64 
    where F: Fn(f64) -> f64 {
        let h = (b - a) / n as f64;
        let mut sum = (f(a) + f(b)) / 2.0;
        
        for i in 1..n {
            let x = a + i as f64 * h;
            sum += f(x);
        }
        
        sum * h
    }
    
    // 辛普森法则
    fn simpson_rule<F>(f: F, a: f64, b: f64, n: usize) -> f64 
    where F: Fn(f64) -> f64 {
        let h = (b - a) / n as f64;
        let mut sum = f(a) + f(b);
        
        for i in 1..n {
            let x = a + i as f64 * h;
            let coefficient = if i % 2 == 0 { 2.0 } else { 4.0 };
            sum += coefficient * f(x);
        }
        
        sum * h / 3.0
    }
    
    // 计算 ∫₀¹ x² dx = 1/3
    let f_x_squared = |x: f64| x * x;
    let exact_integral = 1.0 / 3.0;  // 解析解
    
    println!("数值积分 ∫₀¹ x² dx (解析解 = 1/3 ≈ {:.10}):", exact_integral);
    
    let subdivisions = [10, 100, 1000];
    for &n in &subdivisions {
        let rect_result = rectangle_rule(&f_x_squared, 0.0, 1.0, n);
        let trap_result = trapezoidal_rule(&f_x_squared, 0.0, 1.0, n);
        let simp_result = simpson_rule(&f_x_squared, 0.0, 1.0, n);
        
        println!("  n = {} 个子区间:", n);
        println!("    矩形法则: {:.10}, 误差: {:.2e}", rect_result, (rect_result - exact_integral).abs());
        println!("    梯形法则: {:.10}, 误差: {:.2e}", trap_result, (trap_result - exact_integral).abs());
        println!("    辛普森法则: {:.10}, 误差: {:.2e}", simp_result, (simp_result - exact_integral).abs());
    }
    
    // 计算 ∫₀^(π/2) sin(x) dx = 1
    let f_sin = |x: f64| x.sin();
    let sin_integral_exact = 1.0;
    
    println!("\n数值积分 ∫₀^(π/2) sin(x) dx (解析解 = 1):");
    use std::f64::consts::PI;
    
    let sin_trap_result = trapezoidal_rule(&f_sin, 0.0, PI/2.0, 1000);
    let sin_simp_result = simpson_rule(&f_sin, 0.0, PI/2.0, 1000);
    
    println!("  梯形法则 (n=1000): {:.10}, 误差: {:.2e}", sin_trap_result, (sin_trap_result - sin_integral_exact).abs());
    println!("  辛普森法则 (n=1000): {:.10}, 误差: {:.2e}", sin_simp_result, (sin_simp_result - sin_integral_exact).abs());
    
    // 8.3 泰勒级数逼近
    println!("\n8.3 泰勒级数函数逼近");
    
    // sin(x) 的泰勒展开：sin(x) = x - x³/3! + x⁵/5! - x⁷/7! + ...
    fn taylor_sin(x: f64, terms: usize) -> f64 {
        let mut result = 0.0;
        let mut term = x;  // 第一项是 x
        let mut sign = 1.0;
        
        for n in 0..terms {
            result += sign * term;
            
            // 计算下一项：x^(2n+3) / (2n+3)!
            term *= x * x / ((2 * n + 2) as f64 * (2 * n + 3) as f64);
            sign *= -1.0;
        }
        
        result
    }
    
    // exp(x) 的泰勒展开：e^x = 1 + x + x²/2! + x³/3! + ...
    fn taylor_exp(x: f64, terms: usize) -> f64 {
        let mut result = 1.0;  // 常数项
        let mut term = 1.0;
        
        for n in 1..terms {
            term *= x / n as f64;
            result += term;
        }
        
        result
    }
    
    let test_x: f64 = 0.5;
    println!("泰勒级数逼近 (x = {}):", test_x);
    
    println!("sin(x) 泰勒展开:");
    let sin_exact = test_x.sin();
    for &n in &[3, 5, 10, 15] {
        let sin_approx = taylor_sin(test_x, n);
        let error = (sin_approx - sin_exact).abs();
        println!("  {} 项: {:.10}, 误差: {:.2e}", n, sin_approx, error);
    }
    println!("  精确值: {:.10}", sin_exact);
    
    println!("\nexp(x) 泰勒展开:");
    let exp_exact = test_x.exp();
    for &n in &[5, 10, 15, 20] {
        let exp_approx = taylor_exp(test_x, n);
        let error = (exp_approx - exp_exact).abs();
        println!("  {} 项: {:.10}, 误差: {:.2e}", n, exp_approx, error);
    }
    println!("  精确值: {:.10}", exp_exact);
    
    // 8.4 迭代算法收敛性
    println!("\n8.4 迭代算法收敛性分析");
    
    // 不动点迭代求解 √a：x_{n+1} = (x_n + a/x_n) / 2
    fn fixed_point_sqrt(a: f64, x0: f64, max_iter: usize) -> Vec<f64> {
        let mut iterations = vec![x0];
        let mut x = x0;
        
        for _ in 0..max_iter {
            x = (x + a / x) / 2.0;
            iterations.push(x);
            
            if iterations.len() >= 2 {
                let prev = iterations[iterations.len() - 2];
                if (x - prev).abs() < 1e-15 {
                    break;
                }
            }
        }
        
        iterations
    }
    
    println!("不动点迭代求 √5:");
    let sqrt5_iterations = fixed_point_sqrt(5.0, 2.0, 10);
    let sqrt5_exact = 5.0_f64.sqrt();
    
    for (i, &x) in sqrt5_iterations.iter().enumerate() {
        let error = (x - sqrt5_exact).abs();
        println!("  迭代 {}: x = {:.15}, 误差: {:.2e}", i, x, error);
    }
    
    // 收敛率分析
    if sqrt5_iterations.len() >= 4 {
        let n = sqrt5_iterations.len();
        let e_n = (sqrt5_iterations[n-1] - sqrt5_exact).abs();
        let e_n_minus_1 = (sqrt5_iterations[n-2] - sqrt5_exact).abs();
        let e_n_minus_2 = (sqrt5_iterations[n-3] - sqrt5_exact).abs();
        
        if e_n_minus_1 > 1e-15 && e_n_minus_2 > 1e-15 {
            let convergence_rate = (e_n / e_n_minus_1).log2() / (e_n_minus_1 / e_n_minus_2).log2();
            println!("  估计收敛率: {:.2}", convergence_rate);
        }
    }
    
    println!("\n数值方法总结:");
    println!("  • 牛顿法：二次收敛，但需要导数");
    println!("  • 数值积分：精度随子区间数增加而提高");
    println!("  • 泰勒级数：在收敛半径内精度随项数增加");
    println!("  • 迭代方法：简单但收敛性需要分析");
    
    println!("\n{}", "=".repeat(80));
}

/**
 * 演示实际应用案例
 * 
 * 重点知识：
 * - 科学计算应用
 * - 工程数学问题
 * - 金融数学建模
 * - 信号处理算法
 */
fn demonstrate_practical_applications() {
    println!("\n9. 实际应用案例");
    println!("   浮点数学函数在实际问题中的综合应用");
    
    // 9.1 物理仿真：简谐振动
    println!("\n9.1 物理仿真：简谐振动分析");
    
    // 简谐振动方程：x(t) = A·cos(ωt + φ)
    #[derive(Debug)]
    struct HarmonicOscillator {
        amplitude: f64,      // 振幅 A
        frequency: f64,      // 角频率 ω (rad/s)
        phase: f64,          // 初始相位 φ (rad)
        damping: f64,        // 阻尼系数 γ
    }
    
    impl HarmonicOscillator {
        fn new(amplitude: f64, frequency: f64, phase: f64, damping: f64) -> Self {
            HarmonicOscillator { amplitude, frequency, phase, damping }
        }
        
        // 无阻尼简谐振动位移
        fn position(&self, t: f64) -> f64 {
            self.amplitude * (self.frequency * t + self.phase).cos()
        }
        
        // 有阻尼振动位移：x(t) = A·e^(-γt)·cos(ωt + φ)
        fn damped_position(&self, t: f64) -> f64 {
            self.amplitude * (-self.damping * t).exp() * (self.frequency * t + self.phase).cos()
        }
        
        // 速度：v(t) = -Aω·sin(ωt + φ)
        fn velocity(&self, t: f64) -> f64 {
            -self.amplitude * self.frequency * (self.frequency * t + self.phase).sin()
        }
        
        // 能量分析
        fn kinetic_energy(&self, t: f64, mass: f64) -> f64 {
            let v = self.velocity(t);
            0.5 * mass * v * v
        }
        
        fn potential_energy(&self, t: f64, spring_constant: f64) -> f64 {
            let x = self.position(t);
            0.5 * spring_constant * x * x
        }
    }
    
    use std::f64::consts::PI;
    
    let oscillator = HarmonicOscillator::new(2.0, 2.0 * PI, PI / 4.0, 0.1);
    let mass = 1.0;  // kg
    let spring_constant = (2.0 * PI).powi(2);  // N/m，使得 ω = √(k/m) = 2π
    
    println!("简谐振动参数:");
    println!("  振幅 A = {} m", oscillator.amplitude);
    println!("  角频率 ω = {} rad/s", oscillator.frequency);
    println!("  初始相位 φ = {} rad = {:.1}°", oscillator.phase, oscillator.phase.to_degrees());
    println!("  阻尼系数 γ = {} s⁻¹", oscillator.damping);
    
    println!("\n振动分析 (时间点):");
    let time_points = [0.0, 0.25, 0.5, 0.75, 1.0];
    
    for t in time_points {
        let x = oscillator.position(t);
        let x_damped = oscillator.damped_position(t);
        let v = oscillator.velocity(t);
        let ke = oscillator.kinetic_energy(t, mass);
        let pe = oscillator.potential_energy(t, spring_constant);
        let total_energy = ke + pe;
        
        println!("  t = {:.2}s:", t);
        println!("    位移: x = {:.4} m (无阻尼), {:.4} m (有阻尼)", x, x_damped);
        println!("    速度: v = {:.4} m/s", v);
        println!("    动能: KE = {:.4} J, 势能: PE = {:.4} J, 总能量: {:.4} J", ke, pe, total_energy);
    }
    
    // 9.2 工程数学：傅里叶分析
    println!("\n9.2 工程数学：离散傅里叶变换 (DFT)");
    
    // 简化的DFT实现（仅展示概念）
    fn simple_dft(signal: &[f64]) -> Vec<(f64, f64)> {
        let n = signal.len();
        let mut result = Vec::new();
        
        for k in 0..n {
            let mut real = 0.0;
            let mut imag = 0.0;
            
            for j in 0..n {
                let angle = -2.0 * PI * (k * j) as f64 / n as f64;
                real += signal[j] * angle.cos();
                imag += signal[j] * angle.sin();
            }
            
            result.push((real, imag));
        }
        
        result
    }
    
    // 生成测试信号：f(t) = sin(2πft) + 0.5*sin(4πft)
    let sample_rate = 32.0;  // Hz
    let _duration = 1.0;      // 秒
    let samples: Vec<f64> = (0..32)
        .map(|i| {
            let t = i as f64 / sample_rate;
            (2.0 * PI * t).sin() + 0.5 * (4.0 * PI * t).sin()
        })
        .collect();
    
    println!("信号分析 (前8个采样点):");
    for (i, &sample) in samples.iter().take(8).enumerate() {
        let t = i as f64 / sample_rate;
        println!("  t = {:.3}s: f(t) = {:.4}", t, sample);
    }
    
    // 计算频谱
    let spectrum = simple_dft(&samples);
    println!("\n频谱分析 (前8个频率分量):");
    for (k, &(real, imag)) in spectrum.iter().take(8).enumerate() {
        let magnitude = (real * real + imag * imag).sqrt();
        let phase = imag.atan2(real);
        let frequency = k as f64 * sample_rate / samples.len() as f64;
        
        println!("  f = {:.1} Hz: 幅值 = {:.3}, 相位 = {:.3} rad", frequency, magnitude, phase);
    }
    
    // 9.3 金融数学：Black-Scholes期权定价
    println!("\n9.3 金融数学：Black-Scholes期权定价模型");
    
    // 累积标准正态分布函数的近似
    fn cumulative_normal(x: f64) -> f64 {
        if x < 0.0 {
            1.0 - cumulative_normal(-x)
        } else {
            // Abramowitz and Stegun approximation
            let t = 1.0 / (1.0 + 0.2316419 * x);
            let d = 0.3989423 * (-x * x / 2.0).exp();
            let prob = d * t * (0.3193815 + t * (-0.3565638 + t * (1.7814779 + t * (-1.8212560 + t * 1.3302744))));
            1.0 - prob
        }
    }
    
    // Black-Scholes欧式看涨期权定价
    fn black_scholes_call(s: f64, k: f64, t: f64, r: f64, sigma: f64) -> f64 {
        let d1 = ((s / k).ln() + (r + 0.5 * sigma * sigma) * t) / (sigma * t.sqrt());
        let d2 = d1 - sigma * t.sqrt();
        
        s * cumulative_normal(d1) - k * (-r * t).exp() * cumulative_normal(d2)
    }
    
    // Black-Scholes欧式看跌期权定价
    fn black_scholes_put(s: f64, k: f64, t: f64, r: f64, sigma: f64) -> f64 {
        let d1 = ((s / k).ln() + (r + 0.5 * sigma * sigma) * t) / (sigma * t.sqrt());
        let d2 = d1 - sigma * t.sqrt();
        
        k * (-r * t).exp() * cumulative_normal(-d2) - s * cumulative_normal(-d1)
    }
    
    // 期权希腊字母
    fn option_delta_call(s: f64, k: f64, t: f64, r: f64, sigma: f64) -> f64 {
        let d1 = ((s / k).ln() + (r + 0.5 * sigma * sigma) * t) / (sigma * t.sqrt());
        cumulative_normal(d1)
    }
    
    fn option_gamma(s: f64, k: f64, t: f64, r: f64, sigma: f64) -> f64 {
        let d1 = ((s / k).ln() + (r + 0.5 * sigma * sigma) * t) / (sigma * t.sqrt());
        let n_d1 = (1.0 / (2.0 * PI).sqrt()) * (-0.5 * d1 * d1).exp();
        n_d1 / (s * sigma * t.sqrt())
    }
    
    let stock_price = 100.0;       // 当前股价
    let strike_price = 105.0;      // 行权价
    let time_to_expiry = 0.25;     // 到期时间（年）
    let risk_free_rate = 0.05;     // 无风险利率
    let volatility = 0.2;          // 波动率
    
    let call_price = black_scholes_call(stock_price, strike_price, time_to_expiry, risk_free_rate, volatility);
    let put_price = black_scholes_put(stock_price, strike_price, time_to_expiry, risk_free_rate, volatility);
    let delta = option_delta_call(stock_price, strike_price, time_to_expiry, risk_free_rate, volatility);
    let gamma = option_gamma(stock_price, strike_price, time_to_expiry, risk_free_rate, volatility);
    
    println!("Black-Scholes期权定价:");
    println!("  标的资产价格: ${:.2}", stock_price);
    println!("  行权价格: ${:.2}", strike_price);
    println!("  到期时间: {:.2} 年", time_to_expiry);
    println!("  无风险利率: {:.1}%", risk_free_rate * 100.0);
    println!("  波动率: {:.1}%", volatility * 100.0);
    println!("结果:");
    println!("  看涨期权价格: ${:.4}", call_price);
    println!("  看跌期权价格: ${:.4}", put_price);
    println!("  Delta (价格敏感性): {:.4}", delta);
    println!("  Gamma (Delta敏感性): {:.6}", gamma);
    
    // 验证看涨-看跌平价关系: C - P = S - K·e^(-rT)
    let put_call_parity = call_price - put_price;
    let theoretical_parity = stock_price - strike_price * (-risk_free_rate * time_to_expiry).exp();
    println!("  看涨-看跌平价验证:");
    println!("    C - P = {:.4}", put_call_parity);
    println!("    S - K·e^(-rT) = {:.4}", theoretical_parity);
    println!("    误差: {:.6}", (put_call_parity - theoretical_parity).abs());
    
    // 9.4 信号处理：数字滤波器
    println!("\n9.4 信号处理：简单低通滤波器");
    
    // 一阶RC低通滤波器的数字实现
    struct LowPassFilter {
        cutoff_freq: f64,    // 截止频率
        sample_rate: f64,    // 采样率
        alpha: f64,          // 滤波系数
        prev_output: f64,    // 前一次输出
    }
    
    impl LowPassFilter {
        fn new(cutoff_freq: f64, sample_rate: f64) -> Self {
            let dt = 1.0 / sample_rate;
            let rc = 1.0 / (2.0 * PI * cutoff_freq);
            let alpha = dt / (rc + dt);
            
            LowPassFilter {
                cutoff_freq,
                sample_rate,
                alpha,
                prev_output: 0.0,
            }
        }
        
        fn filter(&mut self, input: f64) -> f64 {
            self.prev_output = self.alpha * input + (1.0 - self.alpha) * self.prev_output;
            self.prev_output
        }
        
        fn reset(&mut self) {
            self.prev_output = 0.0;
        }
    }
    
    // 生成含噪声的信号
    let mut filter = LowPassFilter::new(5.0, 50.0);  // 5Hz截止，50Hz采样
    let signal_freq = 2.0;  // 信号频率2Hz
    let noise_freq = 20.0;  // 噪声频率20Hz
    
    println!("数字低通滤波器 (截止频率: 5Hz):");
    println!("时间   | 原始信号  | 噪声    | 含噪信号  | 滤波输出");
    println!("-------|-----------|---------|-----------|----------");
    
    for i in 0..10 {
        let t = i as f64 / 50.0;  // 时间
        let signal = 3.0 * (2.0 * PI * signal_freq * t).sin();  // 主信号
        let noise = 0.8 * (2.0 * PI * noise_freq * t).sin();    // 高频噪声
        let noisy_signal = signal + noise;                       // 含噪声信号
        let filtered = filter.filter(noisy_signal);             // 滤波输出
        
        println!("{:6.3} | {:8.4} | {:7.4} | {:8.4} | {:8.4}",
                t, signal, noise, noisy_signal, filtered);
    }
    
    // 9.5 统计分析：正态分布和假设检验
    println!("\n9.5 统计分析：正态分布和假设检验");
    
    // 正态分布相关函数
    fn normal_pdf(x: f64, mean: f64, std_dev: f64) -> f64 {
        let coefficient = 1.0 / (std_dev * (2.0 * PI).sqrt());
        let exponent = -0.5 * ((x - mean) / std_dev).powi(2);
        coefficient * exponent.exp()
    }
    
    // Z检验统计量
    fn z_test_statistic(sample_mean: f64, population_mean: f64, population_std: f64, sample_size: usize) -> f64 {
        (sample_mean - population_mean) / (population_std / (sample_size as f64).sqrt())
    }
    
    // 置信区间计算
    fn confidence_interval(sample_mean: f64, std_error: f64, z_critical: f64) -> (f64, f64) {
        let margin_of_error = z_critical * std_error;
        (sample_mean - margin_of_error, sample_mean + margin_of_error)
    }
    
    let population_mean = 100.0;    // 总体均值
    let population_std = 15.0;      // 总体标准差
    let sample_mean = 105.0;        // 样本均值
    let sample_size = 36;           // 样本大小
    let z_critical_95 = 1.96;       // 95%置信区间的临界值
    
    println!("统计假设检验示例:");
    println!("  H0: μ = {} (原假设)", population_mean);
    println!("  H1: μ ≠ {} (备择假设)", population_mean);
    println!("  样本均值: {}", sample_mean);
    println!("  样本大小: n = {}", sample_size);
    println!("  总体标准差: σ = {}", population_std);
    
    let z_stat = z_test_statistic(sample_mean, population_mean, population_std, sample_size);
    let std_error = population_std / (sample_size as f64).sqrt();
    let (ci_lower, ci_upper) = confidence_interval(sample_mean, std_error, z_critical_95);
    
    println!("结果:");
    println!("  Z统计量: {:.4}", z_stat);
    println!("  标准误差: {:.4}", std_error);
    println!("  95%置信区间: [{:.2}, {:.2}]", ci_lower, ci_upper);
    println!("  P值 (双尾): {:.6}", 2.0 * (1.0 - cumulative_normal(z_stat.abs())));
    
    if z_stat.abs() > z_critical_95 {
        println!("  结论: 拒绝原假设 (α = 0.05)");
    } else {
        println!("  结论: 不能拒绝原假设 (α = 0.05)");
    }
    
    // 正态分布概率密度函数绘制（数值形式）
    println!("\n正态分布 N(100, 15²) 的概率密度:");
    let x_values = [70.0, 85.0, 100.0, 115.0, 130.0];
    for x in x_values {
        let pdf = normal_pdf(x, population_mean, population_std);
        let z_score = (x - population_mean) / population_std;
        println!("  x = {:.0}: f(x) = {:.6}, z = {:.2}", x, pdf, z_score);
    }
    
    // 9.6 数值优化：梯度下降
    println!("\n9.6 数值优化：梯度下降法");
    
    // 目标函数：f(x) = x² - 4x + 4 = (x-2)²，最小值在x=2
    fn objective_function(x: f64) -> f64 {
        x * x - 4.0 * x + 4.0
    }
    
    // 导数：f'(x) = 2x - 4
    fn gradient(x: f64) -> f64 {
        2.0 * x - 4.0
    }
    
    // 梯度下降算法
    fn gradient_descent(start: f64, learning_rate: f64, tolerance: f64, max_iterations: usize) -> Vec<(f64, f64)> {
        let mut x = start;
        let mut history = vec![(x, objective_function(x))];
        
        for _i in 0..max_iterations {
            let grad = gradient(x);
            let new_x = x - learning_rate * grad;
            let new_f = objective_function(new_x);
            
            history.push((new_x, new_f));
            
            if (new_x - x).abs() < tolerance {
                break;
            }
            
            x = new_x;
        }
        
        history
    }
    
    println!("梯度下降优化 f(x) = x² - 4x + 4:");
    let optimization_history = gradient_descent(0.0, 0.1, 1e-6, 100);
    
    println!("迭代过程 (前10步):");
    for (i, &(x, f_x)) in optimization_history.iter().take(10).enumerate() {
        println!("  迭代 {}: x = {:.6}, f(x) = {:.8}", i, x, f_x);
    }
    
    if let Some(&(final_x, final_f)) = optimization_history.last() {
        println!("最终结果:");
        println!("  最优解: x* = {:.8}", final_x);
        println!("  最小值: f(x*) = {:.10}", final_f);
        println!("  理论最优: x = 2, f(2) = 0");
        println!("  收敛步数: {}", optimization_history.len() - 1);
    }
    
    println!("\n实际应用总结:");
    println!("✓ 物理仿真：三角函数用于振动分析");
    println!("✓ 工程数学：指数函数用于信号处理");
    println!("✓ 金融建模：对数和指数用于期权定价");
    println!("✓ 统计分析：正态分布函数用于假设检验");
    println!("✓ 数值优化：导数和梯度用于寻找最优解");
    println!("✓ 信号处理：三角函数用于频域分析");
    
    println!("\n{}", "=".repeat(80));
}
